###################################################################
#### Replication Code for                                      ####
#### The Product and Difference Fallacies for Indirect Effects ####
###################################################################

## Loading the Data ##

load(file = "ProdDiffFallacyRepData.RData")

##################################
## Replication Code for Table 1 ##
## Product Fallacy              ##
##################################
print("# Replication of Table 1 #")

## Row 1: Additive Models

print("Row 1: Additive Models")

mod4 <- summary(lm(Democracy ~  Muslim + SNEW125 + Income, data=fd34))
mod4a <- summary(lm(SNEW125 ~  Muslim + Income, data=fd34))

# Avg Effect of IRT on SR
print("Avg Effect of IRT on SR")
print(mod4a$coef[2,1])
print(mod4a$coef[2,2])

# Avg Effect of SR on Dem
print("Avg Effect of SR on Dem")
print(mod4$coef[3,1])
print(mod4$coef[3,2])

# Indirect Effect 
print("Indirect Effect")
print(mod4a$coef[2,1]*mod4$coef[3,1])

B <- 10000
n <- dim(fd34)[1]
est <- rep(0,B)
set.seed(123)
for(b in 1:B){	
	boot.ind <- sample(1:n,replace=TRUE)
	fdtemp <- fd34[boot.ind,]
	mod4a <- lm(SNEW125 ~  Muslim + Income, data=fdtemp)
	mod4 <- lm(Democracy ~  Muslim + SNEW125 + Income, data=fdtemp)
	est[b] <- mod4a$coef[2]*mod4$coef[3]
	}
print(sd(est))

## High/Low Income Split

k <- quantile(fd34$Income,probs=.8)
w <- mean(fd34$Income < k)

## Row 2: High Income
print("Row 2: High Income")
mod4a0 <- summary(lm(SNEW125 ~  Muslim , data=fd34[fd34$Income>=k,]))
mod40 <- summary(lm(Democracy ~  Muslim + SNEW125 , data=fd34[fd34$Income >=k,]))

# Avg Effect of IRT on SR
print("Avg Effect of IRT on SR")
print(mod4a0$coef[2,1])
print(mod4a0$coef[2,2])

# Avg Effect of SR on Dem
print("Avg Effect of SR on Dem")
print(mod40$coef[3,1])
print(mod40$coef[3,2])

# Indirect Effect 
print("Indirect Effect")
print(mod4a0$coef[2,1]*mod40$coef[3,1])

B <- 10000
fd34hi <- subset(fd34,Income>=k)
n <- dim(fd34hi)[1]
est <- rep(0,B)
set.seed(123)
for(b in 1:B){	
	boot.ind <- sample(1:n,replace=TRUE)
	mod4a0 <- lm(SNEW125 ~  Muslim , data=fd34hi[ boot.ind,])
	mod40 <- lm(Democracy ~  Muslim + SNEW125 , data=fd34hi[ boot.ind,])
	est[b] <- mod4a0$coef[2]*mod40$coef[3]
	}
print(sd(est,na.rm=TRUE))

## Row 3: Low Income
print("Row 3: Low Income")
mod4a1 <- summary(lm(SNEW125 ~  Muslim , data=fd34[ fd34$Income < k,]))
mod41 <- summary(lm(Democracy ~  Muslim + SNEW125 , data=fd34[ fd34$Income < k,]))

# Avg Effect of IRT on SR
print("Avg Effect of IRT on SR")
print(mod4a1$coef[2,1])
print(mod4a1$coef[2,2])

# Avg Effect of SR on Dem
print("Avg Effect of SR on Dem")
print(mod41$coef[3,1])
print(mod41$coef[3,2])

# Indirect Effect 
print("Indirect Effect")
print(mod4a1$coef[2,1]*mod41$coef[3,1])

B <- 10000
fd34low <- subset(fd34,Income<k)
n <- dim(fd34low)[1]
est <- rep(0,B)
set.seed(123)
for(b in 1:B){	
	boot.ind <- sample(1:n,replace=TRUE)
	mod4a1 <- lm(SNEW125 ~  Muslim , data=fd34low[ boot.ind,])
	mod41 <- lm(Democracy ~  Muslim + SNEW125 , data=fd34low[ boot.ind,])
	est[b] <- mod4a1$coef[2]*mod41$coef[3]
	}
print(sd(est))

## Row 4:.20/.80 Weighted Average
print("Row 4:.20/.80 Weighted Average")

# Avg IRT on SR effect
print("Avg IRT on SR effect")
mod4a0 <- summary(lm(SNEW125 ~  Muslim , data=fd34[fd34$Income>=k,]))
mod4a1 <- summary(lm(SNEW125 ~  Muslim , data=fd34[ fd34$Income < k,]))
print((1-w)*mod4a0$coef[2,1] + w*mod4a1$coef[2,1])

B <- 10000
n <- dim(fd34)[1]
est <- rep(0,B)
set.seed(123)
for(b in 1:B){	
	boot.ind <- sample(1:n,replace=TRUE)
	fdtemp <- fd34[boot.ind,]
	w.boot <- mean(fdtemp$Income < k)
	mod4a0 <- lm(SNEW125  ~  Muslim , data=fdtemp[fdtemp$Income >= k,])
	mod4a1 <- lm(SNEW125  ~  Muslim , data=fdtemp[fdtemp$Income < k,])
	est[b] <- (1-w.boot)*mod4a0$coef[2] + w.boot*mod4a1$coef[2]
	}
print(sd(est,na.rm=TRUE))

# Avg SR on Dem effect
print("Avg SR on Dem effect")
mod40 <- summary(lm(Democracy ~  Muslim + SNEW125 , data=fd34[fd34$Income >=k,]))
mod41 <- summary(lm(Democracy ~  Muslim + SNEW125 , data=fd34[ fd34$Income < k,]))
print((1-w)*mod40$coef[3,1] + w*mod41$coef[3,1])


B <- 10000
n <- dim(fd34)[1]
est <- rep(0,B)
set.seed(123)
for(b in 1:B){	
	boot.ind <- sample(1:n,replace=TRUE)
	fdtemp <- fd34[boot.ind,]
	w.boot <- mean(fdtemp$Income < k)
	mod40 <- lm(Democracy ~  Muslim + SNEW125 , data=fdtemp[fdtemp$Income > k,])
	mod41 <- lm(Democracy ~  Muslim + SNEW125 , data=fdtemp[fdtemp$Income <= k,])
	est[b] <- (1-w.boot)*mod40$coef[3] + w.boot*mod41$coef[3]
	}
print(sd(est))

# Avg Indirect Effect
print("Avg Indirect Effect")
mod40 <- summary(lm(Democracy ~  Muslim + SNEW125 , data=fd34[fd34$Income >=k,]))
mod41 <- summary(lm(Democracy ~  Muslim + SNEW125 , data=fd34[ fd34$Income < k,]))
mod4a0 <- summary(lm(SNEW125 ~  Muslim , data=fd34[fd34$Income>=k,]))
mod4a1 <- summary(lm(SNEW125 ~  Muslim , data=fd34[ fd34$Income < k,]))

print((1-w)*mod4a0$coef[2,1]*mod40$coef[3,1] + w*mod4a1$coef[2,1]*mod41$coef[3,1])

# Product of Averages (Not in table, but stated in text)
print("Product of Averages (Not in table, but stated in text)")
print(((1-w)*mod4a0$coef[2,1] + w*mod4a1$coef[2,1])*
((1-w)*mod40$coef[3,1] + w*mod41$coef[3,1]))

# SE for Avg Indirect Effect
print("SE for Avg Indirect Effect")
B <- 10000
n <- dim(fd34)[1]
est <- rep(0,B)
set.seed(123)
for(b in 1:B){	
	boot.ind <- sample(1:n,replace=TRUE)
	fdtemp <- fd34[boot.ind,]
	w.boot <- mean(fdtemp$Income < k)
	mod4a0 <- lm(SNEW125  ~  Muslim , data=fdtemp[fdtemp$Income >= k,])
	mod4a1 <- lm(SNEW125  ~  Muslim , data=fdtemp[fdtemp$Income < k,])
	mod40 <- lm(Democracy ~  Muslim + SNEW125 , data=fdtemp[fdtemp$Income > k,])
	mod41 <- lm(Democracy ~  Muslim + SNEW125 , data=fdtemp[fdtemp$Income <= k,])
	est[b] <- (1-w.boot)*mod4a0$coef[2]*mod40$coef[3] + w.boot*mod4a1$coef[2]*mod41$coef[3]
	}
print(sd(est,na.rm=TRUE))



##################################
## Replication Code for Table 2 ##
## Difference Fallacy           ##
##################################
print("# Replication of Table 2 #")

## Row 1: Additive Models
print("Row 1: Additive Models")
mod3 <- summary(lm(Democracy ~  Muslim + Income, data=fd12))
mod4 <- summary(lm(Democracy ~  Muslim + LTDF90P + Income, data=fd12))

# Total Effect
print("Total Effect")
print(mod3$coef[2,1])
print(mod3$coef[2,2])

# Direct Effect
print("Direct Effect")
print(mod4$coef[2,1])
print(mod4$coef[2,2])

# Indirect Effect
print("Indirect Effect")
print(mod3$coef[2,1] - mod4$coef[2,1])

B <- 10000
n <- dim(fd12)[1]
est <- rep(0,B)
set.seed(123)
for(b in 1:B){	
	boot.ind <- sample(1:n,replace=TRUE)
	fdtemp <- fd12[boot.ind,]
	mod3 <- lm(Democracy ~  Muslim + Income, data=fdtemp)
	mod4 <- lm(Democracy ~  Muslim + LTDF90P + Income, data=fdtemp)
	est[b] <- mod3$coef[2] - mod4$coef[2]
	}
print(sd(est))

## Row 2: Interaction Model
print("Row 2: Interaction Model")
mod3 <- summary(lm(Democracy ~  Muslim + Income, data=fd12))
mod4 <- summary(lm(Democracy ~  Muslim * LTDF90P  + Income, data=fd12))
mod4a <- lm(LTDF90P ~  Muslim + Income, data=fd12)

# Total Effect
print("Total Effect")
print(mod3$coef[2,1])
print(mod3$coef[2,2])

# 1st Direct Effect
print("1st Direct Effect")
print(as.numeric((mod4$coef[2,1] + mod4$coef[5,1]*(mod4a$coef[1]+mod4a$coef[2]+mod4a$coef[3]*mean(fd12$Income)))))

# 1st Indirect Effect
print("1st Indirect Effect")
print(as.numeric(mod3$coef[2,1] - (mod4$coef[2,1] + mod4$coef[5,1]*(mod4a$coef[1]+mod4a$coef[2]+mod4a$coef[3]*mean(fd12$Income)))))

# 2nd Direct Effect
print("2nd Direct Effect")
print(as.numeric((mod4$coef[2,1] + mod4$coef[5,1]*(mod4a$coef[1]+mod4a$coef[3]*mean(fd12$Income)))))

# 2nd Indirect Effect
print("2nd Indirect Effect")
print(as.numeric(mod3$coef[2,1] - (mod4$coef[2,1] + mod4$coef[5,1]*(mod4a$coef[1]+mod4a$coef[3]*mean(fd12$Income)))))

# SEs for Direct and Indirect Effects of Row 2
print("SEs for Direct and Indirect Effects of Row 2")

B <- 10000
n <- dim(fd12)[1]
estD0 <- rep(0,B)
estI0 <- rep(0,B)
estD1 <- rep(0,B)
estI1 <- rep(0,B)
set.seed(123)
for(b in 1:B){	
	boot.ind <- sample(1:n,replace=TRUE)
	fdtemp <- fd12[boot.ind,]
	mod3 <- lm(Democracy ~  Muslim + Income, data=fdtemp)
	mod4 <- lm(Democracy ~  Muslim * LTDF90P  + Income, data=fdtemp)
	mod4a <- lm(LTDF90P ~  Muslim + Income, data=fd12)
	estD0[b] <- (mod4$coef[2] + mod4$coef[5]*(mod4a$coef[1]+mod4a$coef[3]*mean(fdtemp$Income)))
	estI0[b] <- mod3$coef[2] - estD0[b]
	estD1[b] <- (mod4$coef[2] + mod4$coef[5]*(mod4a$coef[1]+mod4a$coef[2]+mod4a	$coef[3]*mean(fdtemp$Income)))
	estI1[b] <- mod3$coef[2] - estD1[b]	
	}

# SE for 1st Direct Effect
print("SE for 1st Direct Effect")
print(sd(estD1))
# SE for 1st Indirect Effect
print("SE for 1st Indirect Effect")
print(sd(estI1))
# SE for 2nd Direct Effect
print("SE for 2nd Direct Effect")
print(sd(estD0))
# SE for 2nd Indirect Effect
print("SE for 2nd Indirect Effect")
print(sd(estI0))


#####################################
### Replication Code for Figure 2 ###
### Avg Indirect on Treated       ###
#####################################

print("# Replication of Figure 2 #")

# Women in Government
mod6 <- lm(Democracy ~  Muslim + WGVTNEW  + Income, data=fd56)
mod6a <- lm(WGVTNEW ~  Muslim + Income, data=fd56)
print("Indirect effect implied by additive models (stated in text)")
print(mod6a$coef[2]*mod6$coef[3])



library(MatchIt)
m.out <- matchit(Muslim ~ Income,data=fd56)
m.dat <- match.data(m.out)
m1.dat <- subset(m.dat,Muslim==1)
print("Avg Democracy Score for Muslim Countries (stated in text)")
print(mean(m1.dat$Democracy))
m0.dat <- subset(m.dat,Muslim==0)


## Figure 2 ##
print("Figure 2 produced in a file named WG.pdf")
pdf("WG.pdf")
plot(m1.dat$Income,m1.dat$WGVTNEW,xlab="Log GDP (Income)",ylab="Percent Women in Government (WG)",main="Women in Government, Income, and Islamic Religious Tradition (IRT)",pch=19,ylim=c(0,28))
points(m0.dat$Income,m0.dat$WGVTNEW)
#abline(a=mod4a$coef[1,1]+mod4a$coef[3,1],b=mod4a$coef[2,1],lwd=2)
abline(a=mod6a$coef[1],b=mod6a$coef[3],lty=2,lwd=2)
temp <- m1.dat
temp$Muslim <- 0
z0 <- predict(mod6a,newdata=temp)
points(m1.dat$Income,z0,pch=2)
legend("topleft",legend=c("IRT = 0","IRT = 1","WG predicted for IRT=1 (as if IRT=0)"),pch=c(1,19,2))
dev.off()

mod6 <- lm(Democracy ~  Muslim * WGVTNEW  * Income, data=m.dat)
mod6a <- lm(WGVTNEW ~  Muslim * Income, data=m.dat)

temp <- m1.dat
temp$Muslim <- 0
z0 <- predict(mod6a,newdata=temp)
temp$Muslim <- 1
temp$WGVTNEW <- z0
y1z0 <- predict(mod6,newdata=temp)
print("Residual standard deviation that will be reduced by the restriction to moderate income countries")
print(sd(residuals(lm(WGVTNEW ~ Income, data=m0.dat))))


## Estimates for Moderate Income Muslim Countries
print("Estimates for Moderate Income Muslim Countries (stated in text)")

m.out <- matchit(Muslim ~ Income,data=fd56[fd56$Income<3.6 & fd56$Income>=2.3,])
m.dat <- match.data(m.out)
m1.dat <- subset(m.dat,Muslim==1)
m0.dat <- subset(m.dat,Muslim==0)

mod6 <- lm(Democracy ~  Muslim * WGVTNEW  * Income, data=m.dat)
mod6a <- lm(WGVTNEW ~  Muslim * Income, data=m.dat)

temp <- m1.dat
temp$Muslim <- 0
z0 <- predict(mod6a,newdata=temp)
temp$Muslim <- 1
temp$WGVTNEW <- z0
y1z0 <- predict(mod6,newdata=temp)
print("Estimated Indirect Effect")
print(mean(m1.dat$Democracy) - mean(y1z0))
print("Residual Stanard Deviation (demonstrating slight reduction in variance)")
print(sd(residuals(lm(WGVTNEW ~ Income, data=m0.dat))))

